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Abstract 

An upwind scheme is presented for solving the 
three-dimensional Euler equations on unstructured 
tetrahedral meshes. Spatial discretization is accom- 
plished by a cell-centered finite-volume formulation us- 
ing flux-difference splitting. Higher-order differences 
are formed by a novel cell reconstruction process which 
results in computational times per cell comparable to 
those of structured codes. The approach yields highly 
resolved solutions in regions of smooth flow while avoid- 
ing oscillations across shocks without explicit limiting. 
Solutions are advanced in time by a 3-stage Runge- 
Kutta time-stepping scheme with convergence acceler- 
ated to steady state by local time stepping and implicit 
residual smoothing. 

Solutions are presented for a range of configura- 
tions in the transonic speed regime to demonstrate code 
accuracy, speed, and robustness. The results include an 
assessment of grid sensitivity and convergence acceler- 
ation by mesh sequencing. 

Introduction 

Solution algorithms for the Euler and Navier- 
Stokes equations on unstructured meshes have evolved 
rapidly in recent years, e.g. Refs. 1-12. Key mo- 
tivations behind these developments include the need 
for more geometric flexibility in constructing quality 
meshes around complex configurations, a random data 
structure to better facilitate adapting the mesh to the 
physics of the flow, and a means to more easily ac- 
commodate moving meshes. Potential applications are 
many: more rapid analysis of full aircraft configura- 
tions, airframe component integration, iterative design, 
and store separation to name a few. 

Most of the unstructured algorithms developed to 
date are based on either the finite element method or 
central differencing with added dissipation. Only in 

t Research Engineer, Senior member AIAA 

$ Research Engineer, Member AIAA 

Copyright ©1991 by the American Institute of Aeronautics and Astronautics, 

Inc. No copyright is asserted in the United States under Title 17, U.S. Code. 

The U.S. Government has a royalty-free license to exercise all rights under 
the copyright claimed herein for Governmental purposes. All other rights are 
reserved by the copyright owner. 


recent years has upwind differencing been investigated 
for unstructured grids 7-12 . Upwind differencing uti- 
lizes the propagation of information within a mesh in 
accordance with the theory of characteristics in con- 
structing type-dependent differencing for components 
of the information traveling in opposite directions in 
a separate and stable manner. While this approach is 
more computationally intensive than central differenc- 
ing, it does offer the advantages of being more robust 
and requiring less user interaction. 

This paper describes a new upwind scheme which is 
analytically equivalent to that described in Ref. 12, but 
requires one-half the computational time. The scheme 
is a tetrahedral cell-centered, finite-volume upwind for- 
mulation using flux-difference splitting. Higher-order 
accuracy is achieved by a fast multidimensional linear 
reconstruction algorithm. The approach yields highly 
resolved solutions in regions of smooth flow while avoid- 
ing oscillations across shocks without explicitly apply- 
ing a limiter. Solutions are advanced in time by a 3- 
stage Runge-Kutta time stepping scheme with conver- 
gence accelerated to steady state by local time stepping 
and implicit residual smoothing. 

Governing Equations 

The fluid motion is governed by the time depen- 
dent Euler equations for an ideal gas which express 
the conservation of mass, momentum, and energy for a 
compressible inviscid nonconducting adiabatic fluid in 
the absence of external forces. The equations are given 
below in integral form for a bounded domain f! with a 
boundary dfi 

hi I + J L m) ■* ts -° <» 

where 
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The equations are nondimensionalized with a reference 
density poo and a speed of sound a 0 a . Here h x ,h y , and 
h z are the Cartesian components of the exterior surface 
unit normal fi on the boundary 3fl. The Cartesian 
velocity components are it, v, and w in the x, y, and z 
directions, respectively. The term e 0 is the total energy 
per unit volume. With the ideal gas assumption, the 
pressure and total enthalpy can be expressed as 

P=(7-l)(e 0 -ip(u 2 + v 2 + u> 2 )) 

and 

h 0 = — - + i(u 2 + v 2 + tu 2 ) 

-7-1 P 2 

where 7 is the ratio of specific heats and is prescribed 
as 1.4 for air. 

Equation (1) describes a relationship where the 
time rate of change of the state vector Q within the 
domain Q is balanced by the net flux F across the 
boundary surface 9fl. The domain is divided into a 
finite number of tetrahedral cells, and Eq.(l) is applied 
to each cell. The state variables Q are volume-averaged 
values. It can be shown that the difference equations 
at each cell volume are satisfied exactly when evaluated 
at uniform freestream conditions. 

Spatial Discretization 
Flux Splitting 

Flux quantities are computed using Roe’s 13 flux- 
difference-splitting. The flux across each cell face k is 
computed using Roe’s numerical flux formula 

F k = y 2 [F(Q L ) + F(Qjj)— I A I (Qn - Q l )] k 

Here Qi and Q R are the state variables to the left and 
right of the interface k. The matrix A is computed 
from evaluating 

A = SF/dQ 

with Roe-averaged quantities such as: 


h Q = {h OL + h 0R VPR/PL)/( 1 + Vpr/pl) 

a 2 = (7 - l)(h 0 - (u 2 + t > 2 + w 2 )/ 2 ) 

so that 

F(Q fi ) - F(Q l ) = A[Q* - Qi] 

is satisfied exactly. Introducing the diagonalizing ma- 
trices T and T _1 , and the diagonal matrix of eigenval- 
ues A, then | A | is defined as 

| A |= T | A | f- 1 

The term 


| A | (Qfl — Ql) = T I A I T“ 1 AQ 

in Roe’s flux formula can be reduced to three AF flux 
components, each of which is associated with a distinct 
eigenvalue: 

T | A | T -1 AQ AF X | + | AF 4 | + | AF S | 

with 
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u ± h x a 

V ± flyO, 

w ± h z a 
h 0 ± Ua, 

where U = uh x + vhy + wh z and A U = h x Au+hyAv + 
h z Aw. 

For a first-order scheme, the state of the primitive 
variables at each cell face is set to the cell-centered 
averages on either side of the face. 



P = yJpLPR 
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Modified Higher-Order Scheme 

The challenge in constructing an effective higher- 
order scheme is to determine an accurate estimate of 
the left and right states at the cell faces for the flux 
calculation. Barth and Jespersen proposed a multidi- 
mensional linear reconstruction approach 7 which forms 
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the basis for the present scheme. In the cell reconstruc- 
tion approach, higher-order accuracy is achieved by ex- 
panding the cell-centered solution to each cell face with 
a Taylor series: 

q(s, y, z) = q(*c, y c , Z c ) + Vq c • Ar + 0(Ar 2 ) (2) 
where 

q = [p,u,t>, w,p\ T 


The nodal quantities q n are determined in the 
manner described in Ref. 12. Accordingly, estimates of 
the solution are determined at each node by a weighted 
average of the surrounding cell-centered solution quan- 
tities. It is assumed in the nodal averaging procedure 
that the known values of the solution are concentrated 
at the cell centers, and that the contribution to a node 
from the surrounding cells is inversely proportional to 
the distance from each cell centroid to the node: 


This formulation requires that the solution gradient be 
known at the cell centers. 

A new scheme was proposed in Ref. 12 for esti- 
mating the solution gradient. The general approach 
was to: 1) coalesce surrounding cell information to the 
vertices or nodes of the candidate cell, then 2) apply 
the mid-point trapezoidal rule to evaluate the surface 
integral of the gradient theorem 



over the faces of each tetrahedral cell. Here, Vq denotes 
the volume enclosed by the surface ft. 

It is possible to further simplify the method of Ref. 
12 such that Eq.(3) need not be evaluated explicitly. 
The simplification stems from some useful geometrical 
invariant features of triangles and tetrahedra. These 
features are illustrated for an arbitrary tetrahedral cell 
in Fig. 1. Note that a line extending from a cell-vertex 
through the cell-centroid will always intersect the cen- 
troid of the opposing face. Furthermore, the distance 
from the cell-vertex to the cell-centroid is always three- 
fourths of that from the vertex to the opposing face. 
(For a triangle, the comparable ratio of distance is two- 
thirds). By using these invariants along with the fact 
that Ar is the distance between the cell centroid and 
the face centroid, the second term in Eq.(2) can be 
evaluated sis: 


_ . 3q 

Vq c Ar = — Ar 
or 


+ qn 3 + qn a ) qn 4 


4Ar 


Ar (4) 


Thus, Eq.(2) can be approximated for tetrahedral cells 
by the simple formula: 


q/i.j.j q« V* [ Vs(q«i qnj ^ q«j ) qn<] (s) 

where as illustrated in Fig. 1, the subscripts m, itq, 
denote the nodes comprising face /i, 2,3 of cell c and 
n 4 corresponds to the opposite node. This modified 
scheme is analytically equivalent to that in Ref. 12 and 
results in a factor of two reduction in computational 
time of the flow solver. 


q. = (E^)/(Er) 


(6) 
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where 


= {(x c ,i - x n f + (y e<i - y„) 2 + (*„,< - z n ) 2 ]* 

The subscripts n and c, i refer to the node and sur- 
rounding cell-centered values, respectively. Note that 
the reconstruction process utilizes information from all 
of the cells surrounding the candidate cell, thus pro- 
ducing a truly multidimensional higher-order expansion 
in Eq.(5). For boundary nodes, the surrounding face- 
centered boundary conditions and respective distances 
are used in Eq.(6). 

Upwind schemes generally require the use of lim- 
iter functions to obtain smooth higher-order solutions 
around flow discontinuities, which has been the experi- 
ence with the schemes of Refs. 7-11. While the present 
method has yet to be applied in two-dimensions, experi- 
ence with its application in three-dimensions has shown 
that limiting is not required, and that the method 
correctly captures shocks without oscillations. While 
this unexpected result is beneficial, it seems unlikely 
that a high-order scheme could capture discontinuities 
smoothly without some form of artificial dissipation be- 
ing added. It is quite possible that the particular av- 
eraging procedure employed in the present algorithm 
could add dissipation and locally reduce the accuracy 
across discontinuities. With respect to limiting over- 
shoots in the expansion of Eq.(5), it can be reasoned 
that when the averaging procedure of Eq.(6) is applied 
at a node, the resulting q n represents a weighted mean 
value of the surrounding solution, i.e. q n is bounded by 
the extrema of the surrounding solution. Furthermore, 
in three dimensions the summation of Eq.(6) accesses 
an average of 20 to 22 cells for each node, which re- 
sults in a smoothing of errors introduced from the sur- 
rounding solution. Thus, the expansion will have been 
smoothed and bounded by the procedure and should 
not introduce new extrema into the solution. 

The accuracy of the higher-order scheme has not 
been formally determined. It is known that the Tay- 
lor series expansion and the midpoint-trapezoidal rule, 
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which was applied explicitly in Ref. 12 and implicitly 
in the present algorithm, are both second-order accu- 
rate. Thus, the key to accuracy for this method lies in 
the quality of the averaged solution at the nodes. Sev- 
eral averaging procedures were investigated with known 
test functions, both linear and nonlinear, on an arbi- 
trary tetrahedral grid. Among the alternate procedures 
were those based on cell volumes, and on the inverse dis- 
tances from cell-center to node raised to a power. Error 
was assessed by an RMS average of the local errors rel- 
ative to the exact solution. Volume averaging yielded 
the poorest result, while averaging based on the inverse 
distance raised to the first power (Eq.(6)) worked best. 

Time Integration 

A semidiscrete form of the governing equations 
reads 

V^+R<=°, *= 1,2,3,... (6) 

where 

R* = ^ &Si,j 

3=K(i) 

Ri is the residual accrued by summation of the fluxes 
through the four faces k of a tetrahedral cell i. These 
equations are integrated in time using a fully explicit 
m-stage Runge-Kutta time-stepping scheme developed 
by Jameson et al. 14 : 


q(°) 

= Qr 


q(D 

= q! 0) - 

Q^— R( 0 * 

1 Vi ‘ 

|(m — 1) 
'i 

= q! 0) - 

am-iyr-R) 1 

Q,- m) 

= q! 0) - 

a 

a m y. «•, 

Q n + 1 

= Q,- m) 



where the superscript n denotes the time level, and 
the parenthetical superscripts the stage of the Runge- 
Kutta time stepping. The weighting factors c*i to a rn 
are defined as: 


Qfc = 


1 

m — k + 1 ’ 


k — 1, ..., m 


These values of oc^ will give m-order accuracy in time 
for a linear equation. Preliminary calculations were 
made using both a 3-stage and 4-stage scheme. The so- 
lution and convergence characteristics were essentially 
identical. Thus, a 3-stage scheme was used for the cal- 
culations presented in this paper. 


In many cases, time accuracy in the integration is 
not required. For such cases, the solution convergence 
to steady state is accelerated by local time stepping and 
implicit residual smoothing. 


Local Time Stepping 

Local time stepping accelerates convergence by ad- 
vancing the solution at each cell in time at a CFL num- 
ber near the local stability limit. The expression for the 
local time step was derived with the aid of a 2-D sta- 
bility analysis presented in Ref. 9: 


A ti < v 


Vi 

Ai + Bi + Ci 


(8) 


with 

Ai = (| Ui | +o < )5 i (a:) 

Bi = (| ^ | +ai)^ y) 

Ci = (| w 4 | +ai)5< 2) 

where v is the CFL number, Vi is the cell volume, ai is 
the local speed of sound, and sj x \ and are 
the projected areas of cell i in the x, y, and z directions. 
The local time steps were updated every 50 cycles for 
the results presented in this paper. 


Implicit Residual Smoothing 

The maximum time step can be further increased 
by increasing the support of the scheme through im- 
plicit averaging of the residuals 15 with their neighbors. 
The residuals are filtered through a smoothing opera- 
tor (which is essentially the Laplacian operator for a 
uniform grid): 


R, — Rj + eV^IL 

where 

v 2 S7= £ (r--r^) 

y=«(*) 

The summation uses residuals from the neighboring 
cells which share the faces k with cell i. The resulting 
set of equations can easily be solved by using Jacobi 
iteration 

>'=«(*) i=«(») 

A value for s of about 0.5 is suggested in Ref. 16 to 
maintain a strongly diagonally dominant coefficient ma- 
trix. In practice, two Jacobi iterations are adequate to 
give a good approximation of R~ at all cell centers. 
Residual smoothing was performed during every stage 
of the Runge-Kutta time cycle and resulted in a dou- 
bling of the time step. 
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Boundary Conditions 

For the solid boundaries such as the wing and cen- 
terplane, the flow tangency condition is imposed by set- 
ting the velocities on the boundary faces to their cell- 
center values and then subtracting the component nor- 
mal to the solid surface. Density and pressure bound- 
ary conditions are simply set to the cell-centered value. 
A condition of zero mass and energy flux through the 
surface is ensured by setting the left and right states of 
solid boundary faces equal to the boundary conditions 
prior to computing the fluxes with Roe’s approximate 
Riemann solver. This technique only permits a flux of 
the pressure terms of the momentum equations through 
a solid boundary. 

Characteristic boundary conditions are applied to 
the far-field subsonic boundary using the fixed and ex- 
trapolated Riemann invariants corresponding to the in- 
coming and outgoing waves. The incoming Riemann 
invariant is determined from the freestream flow and 
the outgoing invariant is extrapolated from the interior 
domain. The invariants are used to determine the lo- 
cally normal velocity component and speed of sound. 
At an outflow boundary, the two tangential velocity 
components and the entropy are extrapolated from the 
interior, while at an inflow boundary they sire specified 
as having far-field values. These five quantities provide 
a complete definition of the flow in the far field. 

Results 

A range of results are presented in this section to 
show the speed, accuracy and robustness of the flow 
solver. The speed and the accuracy is shown by way of 
grid sensitivity and mesh sequencing studies while its 
robustness is established by flow solutions on two com- 
plex three-dimensional configurations. All the meshes 
for this study were generated using an improved ver- 
sion of the advancing front grid generation program, 
VGRID3D 17 . 

ONERA M0 Wing 

An ONERA M6 wing has been used for the grid 
sensitivity and the mesh sequencing studies. This con- 
figuration has been widely used as a benchmark to 
evaluate performance of newly developed flow solution 
methods. The wing has a leading edge sweep of 30 
degrees, an aspect ratio of 3.8, a taper ratio of 0.56, 
and symmetrical airfoil sections. The wing has a root 
chord of 0.67 and a semispan 6 of 1.0 with a rounded 
tip. The computational domain is bounded by a rect- 
angular box with boundaries at —6.5 < x < 11.0, 
0.0 < y < 2.5, and — 6.5 < z < 6.5. 

Grid Sensitivity. Transonic solutions were com- 
puted on three grids (Fig. 2) at the same condi- 
tions: Mco = 0.84, and a = 3.06°, to make an assess- 


ment of the grid sensitivity. One of the meshes (Mesh 
l) has cells stretched in the spanwise direction where 
gradients are small while the other two meshes have no 
stretching. Mesh size specifications are listed in Table 
1 . 



Mesh 1 

Mesh 2 

Mesh 3 

Total Cells 
Boundary Faces 

35008 

4046 

108755 

9858 

231507 

16984 

Total Nodes 
Boundary Nodes 

6910 

2025 

20412 

4931 

42410 

8494 


Table 1. Mesh size specifications. 


The computations were performed using the 3- 
stage Runge-Kutta time stepping scheme with local 
time stepping, implicit residual smoothing, and a CFL 
number of 4.0. The solutions were started from 
freestream initial conditions with the first-order scheme 
and run until the L 2 -norm (RMS average of all resid- 
uals) decreased one order of magnitude, at which time 
the solver automatically switched to the higher-order 
scheme. The solution history is plotted in Fig. 3 
against CRAY-2S CPU time to provide a relative com- 
parison of computational effort. Figure 3(a) shows the 
L 2 -norm with a decrease of approximately 2.5 orders of 
magnitude. The convergence history of the lift coeffi- 
cient is shown in Fig. 3(b). Additional details of the so- 
lution characteristics are provided in Table 2. For com- 
parison, computations presented in Ref. 12 for Meshes 
2 and 3 used approximately 3 hours and 8 hours of 
CRAY-2S CPU time, respectively. 



Dimensioned 
Memory, mw 

CRAY-2S 
Time, min. 

Number 

Cycles 

Mesh 1 

2.6 

14 

800 

Mesh 2 

8.1 

89 

1565 

Mesh 3 

17.1 

203 

1716 


Table 2. Solution characteristics. 


A comparison of wing surface pressure contours 
for the three meshes is presented in Fig. 4 with con- 
tour intervals of A(p/p 00 ) = 0.02. Before plotting, the 
computed face-centered boundary quantities were aver- 
aged to the boundary nodes using Eq.(6). The contour 
results show very little overall sensitivity to mesh size. 
As expected, the primary effect of the grid occurs with 
the spatial resolution of the shocks. 

Figure 5 shows the effect of mesh size on the 
streamwise surface C p distribution at six span stations. 
The present results are plotted in comparison to exper- 
imental data at a Reynolds number of 11.7 million 18 . 
The computations on all the three meshes agree well 
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with experiment and demonstrate the accuracy of the 
calculations. The primary effect of mesh size is con- 
fined to regions of large gradients such as the leading- 
edge suction peak and the shock, where the finer mesh 
yields sharper shock definition. Comparisons shown in 
Ref. 12 demonstrated that the solutions for Meshes 
2 and 3 are comparable to those obtained with struc- 
tured codes. Although the deterioration of the solution 
for Mesh 1 is greater at the root and tip stations, it 
should be noted that its solution was obtained with an 
order of magnitude less CPU time and over 6 times less 
memory than that for Mesh 3. 

The force and moment coefficients listed in Table 3 
were computed by integrating the face-centered bound- 
ary pressures. The coefficients for lift, drag, pitching 
moment, and wing root bending moment are based 
on reference quantities of S re f = .5255, l = .67, and 
b re f = 1.0. The pitching moment is referenced about 
the wing apex. 



Mesh 1 

Mesh 2 

Mesh 3 

C L 

.2816 

.2904 

.2911 

C D 

.0141 

.0132 

.0123 

C m 

-.1688 

-.1724 

-.1726 

Crbm 

.1270 

.1283 

.1285 


Table 3. Force and moment coefficients. 


Mesh Sequencing. A solution for a given configu- 
ration can be converged to steady state more efficiently 
by resorting to mesh sequencing. This procedure is a 
convergence acceleration technique in which the con- 
verged flow solution obtained on a coarse mesh is used 
as a starting solution for a finer mesh instead of starting 
from the uniform flow. The final steady state solution 
can be obtained in less overall computational time com- 
pared to one started from freestream initial conditions. 
Of course, interpolation of data between two completely 
different unstructured three-dimensional grids is not 
straight forward. A program which makes use of effi- 
cient octree data structures was written to accomplish 
this task. 

Mesh sequencing was applied to the M6 wing to 
quantify the benefits in terms of reduced computational 
time. The converged solution from the stretched mesh 
(Mesh 1), was interpolated onto the fine mesh (Mesh 3) 
as a starting solution. Figure 6 shows the convergence 
history of the L 2 -norm. To achieve the same level of 
residual reduction of approximately 2.5 orders of mag- 
nitude, the mesh sequenced solution results in a savings 
of 32-percent in CPU time. Even with accounting for 
the interpolation time, the net decrease in CPU time is 
still 29-percent. The final results are identical with and 
without mesh sequencing and, thus, are not shown. 


Low- Wing Transport 

Another quantitative assessment of the flow solver 
is made using a low-wing transport configuration de- 
scribed in Ref. 19. The l/17th-scale configuration con- 
tains a supercritical airfoil and a flow-through represen- 
tation of an advanced turbofan nacelle with a bypass ra- 
tio of approximately 6. The experimental pressure mea- 
surements were obtained in the NASA Langley 16-Foot 
Transonic tunnel 20 at transonic Mach numbers with 
Reynolds numbers in the range of 2.5 x 10 6 based on the 
mean aerodynamic chord of the wing. The present cal- 
culations were made for the condition of M c 0 = 0.768 
and a = 1.116°. 

The computational grid consists of 418,939 cells 
and 75,470 total nodes representing the semispan con- 
figuration. The surface grid (Fig. 7) contains 21,428 
boundary faces and 10,716 boundary nodes, including 
the outer boundaries and plane of symmetry. Grid 
stretching was not applied for these preliminary cal- 
culations. A sufficient definition of the internal flow- 
through nacelle geometry was not available, so the grid 
was terminated at a plane inside the inlet, and at the 
two bypass exit planes. Freestream conditions were pre- 
scribed on the exit plane boundaries. A condition of 
M = 0.Q32M oo was imposed on the inlet plane to bal- 
ance the mass flux. 

Surface pressure contours are presented in Fig. 8. 
Good resolution of the wing shock can be observed 
along with evidence of an inboard lambda shock. 

A comparison of the streamwise C p distributions 
are shown in Fig. 9 for six span stations. The invis- 
cid results are compared with experimental data at a 
slightly higher angle of attack. In general, the agree- 
ment is good and consistent with the expected effects 
of viscosity. The shock is more aft, and the effects of 
flow separation downstream of the shock and in the 
lower-surface cusp region result in a lower experimen- 
tal A C p over the aft region. (It should be noted that 
the stations within 0.463 < r\ < .70 had only three cells 
defining the region between the shock and the trail- 
ing edge.) Similar comparisons are presented in Ref. 
21 for a typical transport configuration with a super- 
critical wing. Comparisons shown between viscous and 
inviscid calculations and transonic experimental data 
show viscous effects comparable in magnitude to those 
observed in the present calculations. 

The solution was obtained using mesh sequencing 
and required approximately 6 hours of CPU time and 
31 megawords of memory. The present grid has a con- 
siderable number of cells clustered in directions of small 
gradients where they are not needed. The total num- 
ber of cells could be greatly reduced with grid stretch- 
ing, which would significantly decrease the amount of 
memory and CPU time required to obtain a satisfac- 
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tory solution. Generation of stretched grids for complex 
configurations is presently an active area of research. 

Space Transportation System (STS) 

A computation was made on the Space Trans- 
portation System (STS) at = 1.05, and 

a — —3.1° to demonstrate the robustness of the flow 
solver in obtaining a solution on a complex geome- 
try with a complex flow field. The semispan grid, 
which includes the orbiter, external tank, and solid 
rocket boosters, consists of 108,538 cells and 21,562 to- 
tal nodes. The surface grid shown in Fig. 10 is repre- 
sented by 13,552 triangular faces and 6,780 nodes, in- 
cluding the outer computational boundaries and plane 
of symmetry. The computations were made with zero 
elevon deflection. 

The solution was obtained with 3410 first-order cy- 
cles with a CFL number of 0.5, then 1275 additional 
cycles at higher-order with a CFL number of 1.0. The 
solution required 235 minutes of CPU time on a CRAY- 
28 and used 8.2 megawords of memory. 

Figure 10 shows a composite picture of the surface 
triangulation and the corresponding pressure contours 
on the full configuration. The centerplane grid is shown 
in Fig. 11 along with the pressure contours in Fig. 
12. The basic features of the flow (shocks, expansions 
etc.) have been well captured in the solution consid- 
ering that only 3-4 layers of fine cells have been used 
close to the body. The quality of the present results, 
which were computed on a relatively coarse grid with 
a relatively small amount of computer time, serves as 
a good demonstration of the flow solver capabilities. 

Code Efficiency 

The preceding solutions were computed on a sin- 
gle processor of the NASA Langley Research Center 
Voyager CRAY-2S. All coding within the flow solver 
portion of the code vectorizes with the standard FOR- 
TRAN compiler. The code statistics are summarized 
in Table 4. 


Dimensioned Memory 
CPU Time, First Order 
CPU Time, Higher Order 


74 words/cell 
18 jus/cell/cycle 
34 //s/cell/cycle 


Table 4. Code statistics. 


For comparison, the prior higher-order algorithm of 
Ref. 12 required 65 ,us/cell /cycle. 

To put these statistics in perspective, structured 
Euler codes generally require from 40 to 50 words of 
memory per cell, and 25 to 35 microseconds of CRAY- 
28 time per cell per cycle for higher-order solutions. 
Thus, the new algorithm yields efficiencies which are 


comparable to those of structured codes. 

To avoid possible confusion, it should be noted 
that a 3-D structured mesh of hexahedral elements con- 
tains the same number of nodes as cells (asymptoti- 
cally), whereas an unstructured mesh of tetrahedral el- 
ements generally contains between 5 and 6 times more 
cells them nodes. As shown in Ref. 12, comparable 
accuracy can be achieved between structured and cell- 
centered unstructured codes with the number of cells 
being of the same order of magnitude. Thus, it is im- 
portant to make comparisons based on the number of 
unknowns computed, i.e. the number of cells for a cell- 
centered scheme. 

Concluding Remarks 

An upwind scheme for solving the three- 
dimensional Euler equations on unstructured tetrahe- 
dral meshes has been presented. The algorithm con- 
sists of a time-explicit cell-centered finite-volume for- 
mulation using flux-difference splitting. Higher-order 
accuracy is achieved by a fast multidimensional linear 
reconstruction algorithm which produces a computa- 
tional efficiency comparable to that of structured algo- 
rithms. The approach yields highly resolved solutions 
in regions of smooth flow while avoiding oscillations 
across shocks without explicitly applying a limiter. 

Results have been presented for a range of config- 
urations at transonic speeds to demonstrate the speed, 
accuracy, and robustness of the flow solver. The poten- 
tial for computational efficiency and accuracy has been 
illustrated by obtaining reasonably good solutions on 
the ONERA M6 wing in 14 minutes of CRAY-2S run 
time with 2.6 megawords of memory. A quantitative 
assessment has been presented for a low-wing trans- 
port configuration to demonstrated the robustness of 
the flow solver in providing an accurate solution on a 
complex geometry. A computation has been made on 
the Space Transportation System at a transonic Mach 
number to demonstrate the robustness of the flow solver 
in obtaining a solution on a very complex geometry 
with a complex flow field. 
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Figure 1.- Geometrically invariant features of 
tetrahedra. 
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(a) Mesh 2. 



(b) Mesh 3. 


Figure 2.- Upper surface mesh for ONERA M6 wing. 
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Mesh 3 
Mesh 1 &; 3 




(a) L 2 -norm (b) Lift coefficient 

Figure 6.- Effect of mesh sequencing on convergence history of ONERA M6 wing. 

M TC = 0.84, a = 3.06°. 










Figure 10.- Surface grid and pressure contours on STS. 
Moo = 1.05, a = -3°. 
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